GitHub

Direct sum

A direct sum is useful when a single unknown consists of several blocks with different meanings, for example

  • a tensor and a scalar,
  • a symmetric tensor and a vector,
  • several tensors with different symmetries.

In such cases, it is often convenient to treat the whole collection as a single object while still keeping the block structure explicit. In Tensorial.jl, this is represented by DirectSumArray. A direct-sum value is constructed with pack, or equivalently with ⊕.

Basic usage

Direct sums are particularly convenient when a problem naturally has several coupled unknowns of different types.

The basic constructors are pack and its alias ⊕.

julia> A = @Mat[1.0 2.0; 3.0 4.0]2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 1.0  2.0
 3.0  4.0
julia> x = pack(A, 3.0)2-element DirectSumVector with storage Float64:
 Space(2, 2)
 Space()
julia> y = A ⊕ 3.02-element DirectSumVector with storage Float64:
 Space(2, 2)
 Space()
julia> x == ytrue

The stored blocks are recovered with unpack:

julia> unpack(x)([1.0 2.0; 3.0 4.0], 3.0)
julia> unpack(x, 1)2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 1.0  2.0
 3.0  4.0
julia> unpack(x, 2)3.0

For symmetric tensor blocks, the internal storage uses Mandel coordinates. This is visible through flatview:

julia> As = symmetric(A)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  2.5
 2.5  4.0
julia> z = As ⊕ 3.02-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> flatview(z)4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 1.0
 3.5355339059327378
 4.0
 3.0

Direct sums and automatic differentiation

One of the main advantages of direct sums is that several coupled variables can be treated as one state while preserving the block structure in derivatives.

As a first example, consider a function of a symmetric tensor block A and a scalar block s:

julia> x = symmetric(@Mat[1.0 2.0; 3.0 4.0]) ⊕ 2.02-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> function f(z)
           A, s = unpack(z)
           dot(A, A) + s * tr(A) + s^2
       endf (generic function with 1 method)
julia> f(x)43.5

The gradient is again a direct-sum object with the same blocks:

julia> g = gradient(f, x)2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> unpack(g)([4.0 5.0; 5.0 10.0], 9.0)

The Hessian is returned as a block matrix:

julia> H = hessian(f, x)2×2 DirectSumMatrix with storage Float64:
 Space(Symmetry(2, 2), Symmetry(2, 2))  Space(Symmetry(2, 2),)
 Space(Symmetry(2, 2),)                 Space()
julia> unpack(H, 1, 1)2×2×2×2 SymmetricFourthOrderTensor{2, Float64, 9}:
[:, :, 1, 1] =
 2.0  0.0
 0.0  0.0

[:, :, 2, 1] =
 0.0  1.0
 1.0  0.0

[:, :, 1, 2] =
 0.0  1.0
 1.0  0.0

[:, :, 2, 2] =
 0.0  0.0
 0.0  2.0
julia> unpack(H, 1, 2)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  0.0
 0.0  1.0
julia> unpack(H, 2, 1)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  0.0
 0.0  1.0
julia> unpack(H, 2, 2)2.0

Thus the derivative already reflects the natural block structure of the problem:

  • unpack(H, 1, 1) is the tensor–tensor block,
  • unpack(H, 1, 2) and unpack(H, 2, 1) are the mixed couplings,
  • unpack(H, 2, 2) is the scalar–scalar block.

The same idea applies to vector-valued maps.

julia> x = symmetric(@Mat[1.0 2.0; 3.0 4.0]) ⊕ 2.02-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> C = symmetric(@Mat[2.0 1.0; 1.0 3.0])2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 2.0  1.0
 1.0  3.0
julia> function F(z)
           A, s = unpack(z)
           B = A + s * C
           t = tr(A) + s^2
           B ⊕ t
       endF (generic function with 1 method)
julia> y = F(x)2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
julia> unpack(y)([5.0 4.5; 4.5 10.0], 9.0)

Its Jacobian is again block-structured:

julia> J = gradient(F, x)2×2 DirectSumMatrix with storage Float64:
 Space(Symmetry(2, 2), Symmetry(2, 2))  Space(Symmetry(2, 2),)
 Space(Symmetry(2, 2),)                 Space()
julia> unpack(J, 1, 1)2×2×2×2 SymmetricFourthOrderTensor{2, Float64, 9}:
[:, :, 1, 1] =
 1.0  0.0
 0.0  0.0

[:, :, 2, 1] =
 0.0  0.5
 0.5  0.0

[:, :, 1, 2] =
 0.0  0.5
 0.5  0.0

[:, :, 2, 2] =
 0.0  0.0
 0.0  1.0
julia> unpack(J, 1, 2)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 2.0  1.0
 1.0  3.0
julia> unpack(J, 2, 1)2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  0.0
 0.0  1.0
julia> unpack(J, 2, 2)4.0

This is precisely the structure needed in many coupled problems.

Linear algebra in direct-sum form

A DirectSumArray also has a flat coordinate representation, which is useful for low-level inspection and for linear algebra operations.

julia> flatview(x)4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 1.0
 3.5355339059327378
 4.0
 2.0
julia> flatview(J)4×4 StaticArraysCore.SMatrix{4, 4, Float64, 16} with indices SOneTo(4)×SOneTo(4):
 1.0  0.0  0.0  2.0
 0.0  1.0  0.0  1.41421
 0.0  0.0  1.0  3.0
 1.0  0.0  1.0  4.0

For symmetric tensor blocks, the flat coordinates use Mandel scaling. This makes the Euclidean inner product of the flat representation consistent with the natural tensor inner product.

In particular, linear solves can be written directly in block form. If J is a DirectSumMatrix and r is a DirectSumVector, one may solve

δx = -J \ r

without manually converting the system to flat coordinates.

Example: return mapping solved by Newton's method

A natural use of direct sums is a local Newton solve in return mapping, where the unknown consists of multiple blocks with different meanings.

As a simple model, consider a small-strain von Mises update written as a nonlinear system for

  • the updated symmetric stress σ, and
  • the plastic multiplier increment Δγ.

We collect them into one direct-sum variable,

julia> σ₀ = SymmetricSecondOrderTensor{3}((1.0, 0.2, 0.1, 0.8, 0.05, 0.5))3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 1.0  0.2   0.1
 0.2  0.8   0.05
 0.1  0.05  0.5
julia> x = σ₀ ⊕ 0.02-element DirectSumVector with storage Float64:
 Space(Symmetry(3, 3),)
 Space()
julia> unpack(x)([1.0 0.2 0.1; 0.2 0.8 0.05; 0.1 0.05 0.5], 0.0)

Suppose we want to solve the local system

\[R(\bm{\sigma}, \Delta\gamma) = \begin{Bmatrix} \bm{\sigma} - \bm{\sigma}^{\mathrm{tr}} + \Delta\gamma \bm{n} \\ \|\operatorname{dev}(\bm{\sigma})\| - (\sigma_y + H\,\Delta\gamma) \end{Bmatrix} = \bm{0},\]

where σᵗʳ is the trial stress, σy is the yield stress, and H is the hardening modulus. The flow direction is given by the associative flow rule,

\[n = \frac{\partial f}{\partial \bm{\sigma}},\]

with f the yield function.

The residual can be written directly in terms of a packed state:

julia> σᵗʳ = SymmetricSecondOrderTensor{3}((2.0, 0.4, 0.2, 1.2, 0.1, 0.9))3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 2.0  0.4  0.2
 0.4  1.2  0.1
 0.2  0.1  0.9
julia> σy = 0.60.6
julia> H = 2.02.0
julia> yield_function(σ, Δγ) = norm(dev(σ)) - (σy + H * Δγ)yield_function (generic function with 1 method)
julia> function R(x)
           σ, Δγ = unpack(x)
           f = yield_function(σ, Δγ)
           n = gradient(σ -> yield_function(σ, Δγ), σ)
           R_σ = σ - σᵗʳ + Δγ * n
           R_γ = f
           R_σ ⊕ R_γ
       endR (generic function with 1 method)
julia> x = σᵗʳ ⊕ 0.02-element DirectSumVector with storage Float64:
 Space(Symmetry(3, 3),)
 Space()
julia> unpack(R(x))([0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0], 0.43279555898864464)

Its Jacobian is obtained automatically:

julia> J = gradient(R, x)2×2 DirectSumMatrix with storage Float64:
 Space(Symmetry(3, 3), Symmetry(3, 3))  Space(Symmetry(3, 3),)
 Space(Symmetry(3, 3),)                 Space()
julia> n = gradient(σ -> yield_function(σ, 0.0), σᵗʳ)3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 0.613222   0.387298    0.193649
 0.387298  -0.161374    0.0968246
 0.193649   0.0968246  -0.451848
julia> unpack(J, 1, 1) ≈ one(SymmetricFourthOrderTensor{3})true
julia> unpack(J, 1, 2) ≈ ntrue
julia> unpack(J, 2, 1) ≈ ntrue
julia> unpack(J, 2, 2) == -Htrue

A Newton step can then be written directly in direct-sum form:

julia> r = R(x)2-element DirectSumVector with storage Float64:
 Space(Symmetry(3, 3),)
 Space()
julia> δx = -J \ r2-element DirectSumVector with storage Float64:
 Space(Symmetry(3, 3),)
 Space()
julia> xnew = x + δx2-element DirectSumVector with storage Float64:
 Space(Symmetry(3, 3),)
 Space()
julia> norm(R(xnew))1.2246705451637024e-16

For the present von Mises example, the Newton update reaches the solution in a single step. This reflects the classical radial-return structure of the problem.

For more general return-mapping problems, however, the local residual is genuinely nonlinear, and several Newton iterations are typically required.

APIs

Tensorial.DirectSumArray — Type
DirectSumArray{Axes, T, N, L}

Container for an element of a direct-sum space stored in Mandel-flat coordinates.

Type parameters

  • Axes: block-space layout along each axis
  • T: scalar storage type
  • N: block-array dimension
  • L: total number of flat stored entries

Notes

  • size(A) returns the block size.
  • flatsize(A) returns the size of the flat Mandel storage.
source
Tensorial.pack — Function
pack(args...)

Construct a DirectSumArray from the blocks of a direct-sum space.

Each argument is interpreted as one block of the direct sum. Scalar blocks are stored directly. Tensor blocks are flattened and concatenated into the internal storage. If a block belongs to a symmetric tensor space, Mandel coordinates are used internally.

⊕ (typed by \oplus<tab> ) is provided as an alias for pack.

See also: unpack, flatview

Examples

julia> A = @Mat[1.0 2.0; 3.0 4.0]
2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 1.0  2.0
 3.0  4.0

julia> x = pack(A, 3.0)
2-element DirectSumVector with storage Float64:
 Space(2, 2)
 Space()

julia> As = symmetric(A)
2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  2.5
 2.5  4.0

julia> v = @Vec[4.0, 5.0]
2-element Vec{2, Float64}:
 4.0
 5.0

julia> y = As ⊕ 3.0 ⊕ v
3-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()
 Space(2,)
source
Tensorial.unpack — Function
unpack(A::DirectSumArray)

Return the blocks stored in a DirectSumArray.

unpack(A) returns all direct-sum blocks as a tuple. unpack(A, I...) returns the block at block index I. Tensor blocks are reconstructed from the internal flat storage. For symmetric tensor blocks, the internally stored Mandel coordinates are converted back to the corresponding tensor values.

At the level of block values, unpack is the inverse of pack.

Info

unpack(A) is type-stable, but unpack(A, I...) may be type-unstable if the block index is not known at compile time, since different block indices may correspond to different return types. It is therefore most suitable when I is a constant and constant propagation can resolve the selected block.

See also: pack, flatview

Examples

julia> A = symmetric(@Mat[1.0 2.0; 3.0 4.0])
2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  2.5
 2.5  4.0

julia> x = pack(A, 3.0)
2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()

julia> unpack(x)
([1.0 2.5; 2.5 4.0], 3.0)

julia> unpack(x, 1)
2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  2.5
 2.5  4.0

julia> unpack(x, 2)
3.0
source
Tensorial.flatview — Function
flatview(A::DirectSumArray)

Return A as a flat coordinate array. Symmetric tensor blocks are represented in Mandel coordinates.

See also: pack, unpack

Examples

julia> A = @Mat[1.0 2.0; 3.0 4.0]
2×2 Tensor{Tuple{2, 2}, Float64, 2, 4}:
 1.0  2.0
 3.0  4.0

julia> x = pack(A, 3.0)
2-element DirectSumVector with storage Float64:
 Space(2, 2)
 Space()

julia> flatview(x)
5-element StaticArraysCore.SVector{5, Float64} with indices SOneTo(5):
 1.0
 3.0
 2.0
 4.0
 3.0

julia> As = symmetric(A)
2×2 SymmetricSecondOrderTensor{2, Float64, 3}:
 1.0  2.5
 2.5  4.0

julia> y = pack(As, 3.0)
2-element DirectSumVector with storage Float64:
 Space(Symmetry(2, 2),)
 Space()

julia> flatview(y)
4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 1.0
 3.5355339059327378
 4.0
 3.0
source